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Abstract 

'We introduce a method of estimating parameters associated with a fractal random scattering medium, which utilizes the 
multiscale properties of the scattered field. The example of ray-density fluctuations beyond a phase screen with fractal 
slope is considered. An exact solution to the forward problem, in the case of the Brownian fractal, leads to an expression 
for the volatility of the slope. This expression is invariant under a change of probability measure, a fact which gives 
Tise to the corresponding result for a (stationary) Ornstein-Uhlenbeck slope. We demonstrate that our analytical results 
^Hare consistent with numerical simulations. Finally, an application to the determination of sea ice thickness via sonar is 
discussed. 
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.1. Introduction 



When a wave is scattered by a random medium the 
. scattered radiation can be described by a random field. 
^ We consider the problem of estimating parameters related 
• i-H ,to the characteristics of the scattering medium. The nu- 
^^^■nierous applications of such problems, for example, in the 
(— I areas of industry, geophysics and medicine are well known, 
Q^and are of considerable practical importance 

One possible approach is the following: First, derive 
an expression for the "random scattered field", or some 
J> observable properties thereof. Then, compare these theo- 
Tetical quantities with the corresponding sample quantities 
derived from experiment. A popular version of this strat- 
egy is the "method of moments" which uses the mean, 
yariance, correlation function etc. A refinement of this, 
also along Bayesian lines, is the method of maximum like- 
ilihood, which can provide optimal estimates of the pa- 
rameters in question However, this technique is not 
iusually available in scattering problems. Elementary con- 
. ^ siderations soon lead to the conclusion that scattered fields 
[(arising from, for example, a homogeneous random scat- 
tering medium) are rarely endowed with convenient prop- 
erties. For instance, in the example of surface scattering 
the Markov property will not hold in general, since the 
scatter could be a function of the whole of the scattering 
surface. 

If the scattered field has fractal properties, there is an- 
other possibility for estimating parameters, unrelated to 
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Bayes theorem. Information can be derived from fiuctu- 
ations in the field at all scales. The principle can be un- 
derstood by relation to the "well known fact" that the 
volatility of an Ito diffusion can be recovered with proba- 
bility one_;_given a sample path over any non-zero interval, 
see e.g. [3|. It is common for naturally occurring struc- 
tures to possess multiscale character suitable to be mod- 
eled as a fractal. Moreover, it seems reasonable to suppose 
that such structures might confer multiscale properties to 
the scattered field. Parameter estimation based around 
this idea is potentially very efficient, since it does not re- 
quire the availability of data over many correlation lengths. 
This approach to sensing (inverse scattering) is apparently 
new, and is illustrated here through the example of ray- 
density fluctuations beyond a subfractal phase-changing 
screen (SPS). The concept of rays is of great utility in 
the theory of scattering, most famously in the shortwave 
limit [3l , but also more generally Q . The ray-density, de- 
noted by R, describes intensity fluctuations induced by a 
SPS in an incoherent conflguration (no interference). R 
also features in expressions for the moments of the scat- 
tered intensity in a coherent configuration We con- 
sider scatter from a one-dimensional SPS in an incoherent 
configuration, which provides the simplest example of the 
proposed technique. 

2. The phase screen model 

When a plane of parallel rays passes through a non- 
flat refracting layer, differences in the refractive index of 
the layer cause the rays to scatter. The situation may be 
described by specifying the slope (gradient) of a surface 
of constant phase in terms of a one-dimcnsional stochastic 
process {Yt)t£R- If Y is fractal, then we have a "subfrac- 
tal phase-changing screen" Q. This model also describes 
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surface scattering, when the source and point of observa- 
tion are coincident, and shadowing and multiple scattering 
are neglected, see Fig. [TJ These surfaces can be thought 
of as possessing a faceted structure, and have been used 
to model, among other things, scattering of radio waves 
in the ionosphere Q, and the effect of internal waves on 
acoustic propagation in the ocean Q. 
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Figure 1: Rays scattered from a faceted random surface, or through 
a subfractal phase-changing screen. 

The ray-density can be approximated at a G An 
[a„, a„+i) where a„ = nAa, a unit distance from the layer 
by {R^(At, Aa); n e Z} defined as 
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where tk = kAt and Aa ^ At. The ray-density {R^)a£^ 
is defined as the limit (when it exists) 



Ri := hm hm Ri^ (At, Aa) 



dtSiYt-t + a) (2) 



where a € An. Assume that Y is an Ito diffusion described 
by the stochastic differential equation dVt ~ —9Ytdt + 
(jdi?t, where {Bt)teR is a Brownian motion, and 6 G [0, oo) 
and (T > arc parameters known as the damping and 
volatility respectively. Norris Q found an exact solution 
to the forward problem when 9^0, obtaining when 
Yf = aBf. This is referred to as the case of "Brownian 
slope" . 

3. Brownian slope 

The result from Q, that is an Ito diffusion when 
Yt = aBt is stated briefiy below: 
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known as the occupation density or local time, measures 
infinitesimally the "time spent at x by X before f . Set 
Xt = aBf — t, where {Bt)t>o is a Brownian motion, started 
at zero. Also, let T denote the first time that X hits 
level K < 0. It follows from the Ray-Knight theorem 
on Brownian local time, that the process {L^)x>k has 
generator 



^(^)%(l(.,o)-0 



(4) 



Since Xt — —oo as t — > 
Fx := limt_j.oo < co 
a diffusion, and i?^^ ~ 
a G M; with generator 



oo, {Xt)t>o has a final local time 
. Taking K -cx), {Fx)xm is 
lim;,^oo F-a-b is a diffusion on 
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Now consider the problem of finding a. Calling our 
probability measure P, it follows immediately from (|4]) that 
Lt satisfies 



4 dxL-T/{L)j.i 



(6) 



where I > J > K , and {L)jj denotes the quadratic varia- 
tion of Lt on [J, I]. Also, by sending K — > —oo in or 
from (O, R""^ satisfies 



4 1^ daRa/{R)i,.j 



(7) 



for I < J. At first glance, the expression ([7]) might appear 
to be the end of the matter, but there is an important 
sense in which it is unphysical. A density is of course 
an idealization. In practice, one seeks an approximation 
by measuring the energy flux over some finite area. For a 
general density {Zx)xm deflne the "discrete local average" 

{Zn)n£Z by 

Zn :— — / dxZx (8) 
Aa Ja„ 

where Aa > is the width over which the average is taken. 
The quantity Z is observable, and we assume that mea- 
surements of Z arc available. Consider the case Z = B, a. 
Brownian motion. One quickly finds that 



[{Bn+i - Bnf] = 2Aa/3. 



Define a partition of [/, J] as {a,; 
quadratic variation of Z as 
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{Z)j,j := lim Qij{Z) 

Aa— s-0 



(9) 

, N}, and the 
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where 
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It follows from © that {B)ij ^ (B) ij. If (B) ij is deter- 
ministic dH) implies that = 2(J - /)/3 P-a.s.. This 
follows from the second moment method if 



lim E[Q|„,(B)] 
= lim E[Q/.j(B)]^ 

Aa— ^0 



4(J- 7)2/9. 



(12) 
(13) 
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The expectation in ([12)) is not as easy to evaluate as in 
the Brownian case (without the local average), since the 
increments -B„_|_i — Bn are not independent. However, the 
(quadruple) integrals that result from multiplying out the 
squared sum are trivial to compute using the correlation 
function E[B^BtB„S„] = s{2t + u) where Q<s<t<u< 
V. For example, for hq ~ 0, N ~ 1 we have 
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E[Q\B)]=E[{B,^Bo)' 



=E [B^ - 2BlB^ 



&BlBl - 2BoBf + Sf] 



(14) 



Each of the terms in (|14p can be evaluated by exchanging 
the order of the expectation and the integral, and time 
ordering. For example, the first term: 



Figure 2; Three values of a computed from II16II using 1^ for Brow- 
nian slope. Note that the rate of convergence is only limited by the 
finiteness of At and Aa. The Brownian motion is simulated exactly 
with Aa = 2 X 10"'' and At = 10"*. 



P must also have probability one w.r.t. Q. Let {Bt)t>o be 
a P-Brownian motion, started at zero. Then {Mt)o<t<T 
where 



■ Mt= (17) 

=(Aa)V3 



(15) 



For larger values of N one can use a computer to sum 
the terms and verify that E[(5f j{B)] does indeed appear 
to tend to the required limit ([T5|. Taking this as an as- 
sumption, and applying the same reasoning to the process 
defined by ([S]), it follows from (O, in the limit Aa — >■ 0, 
that satisfies 



daRa/{R) 
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is a martingale, which can be verified, for example, with 
the help of Example 3, p. 233, [l3|. Now, for the process 
(Xt)o<t<T defined as above, it follows from the Girsanov 
theorem (see for example Q) that 

dXt = -l{t<_K}d{Xt + t)dt -dt + crdWt w.r.t. Q (18) 

where {Wt)o<t<T is a Q-Brownian motion. The process 
{Ut)o<t<-KA,T where Ut := Xt+t is a Q-OU process, such 
that dUt = —OUtdt + adWt w.r.t. Q. By the equivalence 
of P and Q, © implies that 



where := limAa^o = Ra fo^' a £ An. Hence (|T6l) 
solves the inverse problem for the Brownian slope. Noting 
that Rj^ = limAt^o R^(At, Aa), R can be simulated using 
(dl . Fig. [5] shows values of a computed using . 

The Brownian slope model is somewhat unsatisfactory 
since Brownian motion is non-stationary and limt_^oo l^t I = 
oo. Of greater verisimilitude is the model resulting when 
9 > 0, then F is a stationary Ornstein-Uhlenbeck (OU) 
process, a Gaussian Markov process with exponential au- 
tocorrelation function. This is referred to as the case of 
"OU slope" and is considered next. 

4. OU slope 

Define a probability measure Q such that dQ := A/tcIP. 
If {Mt)o<t<T is a martingale then P and Q are equivalent 
which means that any event that has probability one w.r.t. 



dxL^r/{L)jj 



-a.s. 
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where I > J > K. Note that Lt is not Markov under Q. 
The remaining steps are understood to be w.r.t. Q. Taking 
K — —oo, {Xt)t>o has a final local time < oo. Then, 
R^ = limf,_j.oo F-a-b where {Ut)te& is a stationary OU 
process described by dlJt = —OUtdt + adWt, and (Wt)tgR 
is a Brownian motion. It now follows from (jl9p . in this 
limit, that R^ satisfies 



a" =4 



daRal{R)u:j 



-a.s. 



(20) 



for I < J. Then, noting that the change-of-measure ar- 
gument goes through as above, and assuming (jl3p . R^ 
satisfies 



daRa/{R) 



I, J 



(21) 



where Ra is defined as in ((16)) . Hence cr can be obtained 
given the ray-density on any non-zero interval in the case 
of OU slope, see Fig. [3] We emphasize that the expression 
is valid for arbitrary 9 € [0, oo). 
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Figure 3; Three values of a computed from II21I I using 1^ for OU 
slope. As in the Brownian case, the rate of convergence is only 
limited by the finiteness of At and Aa. The OU process is simulated 
exactly with Aa = 2 X 10"" and At = 10~*. 

With regard to 6, consider the moments of . It is 
not hard to see that the moments of R^ are equal to the 
moments of R^ in the limit Aa — > 0. Therefore, for ease 
of calculation, the moments of R^ are considered, and the 
superscript is dropped to aid clarity. One quickly finds 
that E[i?] = 1. For the second moment, define the process 
{Ra{to);ae M} as 



Ra{to) 



dtS{Ut-t + a). 



to 



where U is an OU process with transition density 



p{Ut^y\Us^x) = 

V0 



(22) 
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t-s|^Jl/2 



exp 



0{y — xe 



cr2(i - e-2e|t-s|) 



see for example For Xt := Ut — t, it follows that 
E[Rl{to)\{Xt, ^ xo}] (24) 

/OO 
dqp{Ug+a ^q + to\Uo^ UQ)E[Rg{to)\{Xt, - q}]. 
-a 

The unconditional second moment can be found by taking 
to — > — oo, or, equivalently, by fixing to = and sending 
a — > oo: 

E[i?2] = lim E[Rl{0)\{Xo ^ q}] 
=2 [ dqp{U = q)E[R,{0)\{Xo - q}] (25) 



where p(U = q) is the point-density of the long-term limit 
of U. This calculation is similar to the Kac moment for- 
mula In general (|25p must be evaluated numerically. 
An approximate analytical result is available when the cor- 
relation length of U is small compared to the height of the 
point of observation above the surface, that is 1/9 <C 1. 
Assume for simplicity that cr^ ~ 1, then 



lim E\R^ 

Aa^O 



]E[i?2 



(26) 



holds to a close approximation for 9 > 10^, see Fig. S) Fi- 
nally, (|26p combined with (PT|) provides a way to estimate 
of 9 given a. 




Figure 4: (color online) An approximate expression for the second 
moment. Solid lines are I I26II . The denote exact values from 
numerical integration of II25I I. 



5. Application to sensing of sea ice thickness and 
discussion 

An example of an application is the following: The first 
convincing evidence for the thinning of the Arctic sea ice 
was collected using submarine mounted sonar [l2| . Sonar 
methods continue to this day to provide the most reli- 
able source of information on sea ice thickness on large 
scales 131. The thickness of the ice can be inferred from 



the range of the mean surface of the underside. This sur- 
face is known to possess fractal properties, and to have an 
approximately exponential correlation function [l^ . The 
statistics of the ice surface are important parameters for 
determining the relation between the range and the time 
elapsed before the earliest return [iBj . Typically, returns 
from various reflecting points are resolvable in time. It 
follows that the individual intensities can be added, the 
sum being proportional to the ray-density, assuming the 
geometrical optics approximation. Taking the OU process 
as a one-dimensional model of the slope of the surface. 
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estimates of {(7,6} can be obtained given as outlined 
above. 

To summarize, the example provided exploits fluctua- 
tions in the ray-density, which is proportional to the inten- 
sity in incoherent configurations. The potential scope of 
the method is clearly far wider than the example given. As 
already mentioned, the ray-density features in expressions 
for the moments of the intensity in coherent configura- 
tions, in the short wavelength limit. Furthermore, it is R 
which appears to dictate the structure of intensity fluctu- 
ations on a small-scale when the outer scale is small 
(for OU slope this means that 1/9 1). For finite wave- 
length, the requirement that the scattered amplitude must 
satisfy the Helmholtz equation means that the field cannot 
be "nowhere differentiable" , and will therefore have zero 
quadratic variation. However, it is only necessary for frac- 
tal structure to be present down to the resolution of the 
measurements . 
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